void fft (Forward, Size, VecRe, VecIm)
int   Forward, Size;
float * VecRe, * VecIm;
{
   int    LenHalf, Stage, But, ButHalf, i,j,k, ip, pot2;
   float  ArcRe,ArcIm, dArcRe,dArcIm, ReBuf,ImBuf, ArcBuf;
   double Arc;
   pot2 = 0;
   while (Size != (1 << pot2))  pot2++;
   LenHalf = Size >> 1 ;
   j = 1;
   for (i=1; i<Size; i++) {
      if (i<j) {
         ReBuf = VecRe[j-1];
         ImBuf = VecIm[j-1];
         VecRe[j-1] = VecRe[i-1];
         VecIm[j-1] = VecIm[i-1];
         VecRe[i-1] = ReBuf;
         VecIm[i-1] = ImBuf;
      }
      k = LenHalf;
      while (k<j) {
         j -= k;  k = k >> 1;
      }
      j += k;
   }
   for (Stage=1; Stage<=pot2; Stage++) {
      But = 1 << Stage;
      ButHalf = But >> 1;
      ArcRe = (float)1;
      ArcIm = (float)0;
      Arc = (double) (PI/ButHalf);
      dArcRe = (float) cos(Arc);
      dArcIm = (float) sin(Arc);
      if (Forward)  dArcIm = -dArcIm;
      for (j=1; j<=ButHalf; j++) {
         i = j;
         while (i<=Size) {
            ip = i + ButHalf;
            ReBuf = VecRe[ip-1] * ArcRe - VecIm[ip-1] * ArcIm;
            ImBuf = VecRe[ip-1] * ArcIm + VecIm[ip-1] * ArcRe;
            VecRe[ip-1] = VecRe[i-1] - ReBuf;
            VecIm[ip-1] = VecIm[i-1] - ImBuf;
            VecRe[i-1]  = VecRe[i-1] + ReBuf;
            VecIm[i-1]  = VecIm[i-1] + ImBuf;
            i += But ;
         }
         ArcBuf = ArcRe;
         ArcRe = ArcRe  * dArcRe - ArcIm * dArcIm;
         ArcIm = ArcBuf * dArcIm + ArcIm * dArcRe;
   }  }
   if (Forward) {
      for (j=1; j<=Size; j++) {
         VecRe[j-1] /= Size;
         VecIm[j-1] /= Size;
}  }  }
